Effective shape and phase behaviour of short charged rods 
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We explicitly calculate the orientation-dependent second virial coefficient of short charged rods 
in an electrolytic solvent, assuming the rod-rod interactions to be a pairwise sum of hard-core and 
segmental screened-Coulomb repulsions. From the parallel and isotropically averaged second virial 
coefficient, we calculate the effective length and diameter of the rods, for charges and screening 
lengths that vary over several orders of magnitude. Using these effective dimensions, we determine 
the phase diagram, where we distinguish a low-charge and strong-screening regime with a liquid 
crystalline nematic and smectic phase, and a high-charge and weak-screening regime with a plastic 
crystal phase in the phase diagram. 



I. INTRODUCTION 

The study of suspensions of non-spherical colloidal par- 
ticles started with the experimental works of Zocher [l| 
and Bawden et al. Q, and with Onsager's theoretical 
work 3]. It has since developed into a very versatile 
field of research. A lot of attention has been focussed on 
needle-shaped rods, either naturally occurring ones such 
as viruses like Tobacco Mosaic Virus or fd virus 0, BL 3] , 
or laboratory-synthesized ones such as Boehmite rods 6| . 
In recent years, however, a plethora of non-spherical par- 
ticles have been synthesized that are not extremely elon- 
gated, for example: ellipsoidal colloids with aspect ratio 
~ 3 0; colloidal dumbbells 8]; or nano-particles with 
the shape of a rod, disk, snowman, cube, cap, or rasp- 
berry, |t| EE El EE El El El- These particles are 
often charged when dissolved in a polar solvent like wa- 
ter, and hence their pair-interactions involve not only the 
anisotropic steric short-range repulsions but also electro- 
static long-range repulsions. The strength of the latter is 
determined by the charge on the particle and the range 
is determined by the Debye screening length of the sol- 
vent [HI [TtJ • For small charges and strong screening 
(i.e. high salt concentrations), one expects the steric in- 
teractions to be dominant (if we assume that dispersion 
forces can be neglected). Hence, one can use computer 
simulations or theoretical studies of hard anisotropic bod- 
ies [H, [H HE HH, HE to get an idea of the phase 
diagram of the system as a function of concentration. 
In the case of a high charge or weak screening (i.e. low 
salt concentration), however, the situation is less clear- 
cut. There, the degree of anisotropy of the electrostatic 
interactions is not obvious from the outset: on the one 
hand one expects the soft screened-Coulomb interactions 
to wash out the hard-core anisotropy such that the inter- 
actions become effectively more spherically symmetric, 
while on the other hand there are the intriguing findings 
reported for example in Refs. [H HH- The studies in 
these papers apply to systems of charged anisotropic par- 
ticles in a screening medium. It was found that the elec- 
trostatic anisotropy persists to infinitely large distances 



as the asymptotic decay of each multipole contribution to 
the electrostatic potential due to a nonspherical charge 
distribution is equal [HI, HH . This conclusion is in sharp 
contrast to the case of a charge distribution in vacuum, 
where the monopole potential decays more slowly than 
that of the dipole, as each order multipole contribution 
decays slower than the next one does. In our paper we 
investigate the interplay between hard-core and electro- 
static interactions for non-spherical particles, for the rel- 
atively simple particle shape of spherocylinders. 

It is well-established by now that non-spherical col- 
loidal particles can form a wealth of phases in thermo- 
dynamic equilibrium. Nccdlc-likc colloidal rods, for in- 
stance, form a phase sequence I-N-Sm-X upon increas- 
ing the concentration from very dilute up to close pack- 
ing, where I is the completely disordered isotropic fluid 
phase, N the liquid crystalline nematic phase with ori- 
entational ordering, Sm the smectic-A phase built from 
orientationally ordered liquid-like la yers , and X a fully or- 
dered crystal phase 0, HTl EE HE 113, HE HE HJ . This 
phase sequence for colloidal needles is well-established 
for hard-core interactions [H, HE H3]- Also, for softer 
electrostatic screened-Coulomb repulsions in the case of 
charged needles, at least in the regime where the length of 
the rods is much larger than the diameter and the screen- 
ing length of the electrolytic solvent. This ensures that 
the effective diameter of the rods is much smaller than the 
length [J,l3l|. By contrast, particles with shapes that are 
sufficiently close to spherical are not expected to exhibit 
the liquid crystalline phases N and Sm due to their small 
anisotropy. Instead, for such near-spherical particles one 
would expect a plastic crystalline phase (P) to appear 
in the phase diagram, residing in between the isotropic 
fluid and the fully ordered crystal. The P phase is char- 
acterized by positional ordering on a lattice, but with- 
out long-ranged orientational ordering of the particles. 
For instance, a phase sequence I-P-X upon increasing 
the concentration has indeed been established in simula- 
tions of short hard spherocylinders and of hard dumbbells 
with a leng th-to-diameter ratio smaller than about 0.35 
[Til . HE |33| | The question we address in this paper con- 
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cerns the effect of colloidal charge and ionic screening on 
the effective shape of relatively short rods, and on their 
expected phase sequence upon increasing the concentra- 
tion. On the basis of the well-established increase of the 
effective diameter of charged needles compared to their 
hard-core diameter [3l| , it is to be expected that high col- 
loidal charges and weak-screening conditions (i.e. low salt 
concentrations) lead to a decreased anisotropy of short 
charged rods. Hence, this will lead to a larger tendency 
of the system to exhibit a plastic crystal phase instead 
of liquid crystalline phases in the phase diagram, even 
if the hard-core shape would allow for liquid crystalline 
equilibrium phases. 

Of course, suspensions of charged rods have been ex- 
tensively studied theoretically before. Many of these 
studies are based on Onsager's second virial theory for 
hard rods [H, which is modified and extended to take 
into account the effects of charge and screenin g on the 
isotropic-to-nematic transition [3l|, [34], [H, l36l. I37L l38l . 
Some of these studies, for example those of Refs. [3, [3l], 
34], focus on the needle limit in which the rod length 
is very large compared to the screening length. In this 
limit, only the diameter is affected by the electrostatic 
effects, but in such a way that the effective geometry of 
the rod remains needle-like. In Refs. [3a, |36[ rod lengths 
of the order of (or larger than) the Debye length are 
considered, at the expense, however, of ignoring many 
of the prefactors such that the theory is essentially a 
scaling theory. Interestingly, this scaling theory predicts 
nematic-nematic coexistence in some parameter regime, 
which was later confirmed in Ref. (37[. This coexistence 
regime is characterized by a small rod charge density, 
such that the effective geometry of the rod is no longer 
necdlc-likc. Another limit that was studied in detail is 
the limit of weak electrostatic interactions, which nat- 
urally leads to a perturbative description [53, [M, EJ- 
These schemes are very successful at describing the ef- 
fective (non-needle-like) geometry that shows up in the 
angular dependence of the second virial coefficient. An- 
other very interesting effect was identified in Ref. [38[, 
where the correlation free energy of the many-body sys- 
tem of charged rods and counterions was calculated, re- 
sulting in an enhanced tendency to orientational ordering 
and also the possibility of nematic-nematic coexistence. 
With the notable exception of Ref. [4(| , however, most of 
these works on charged rods focus on the isotropic and 
nematic phases and hence, implicitly, on rods which are 
sufficiently elongated to give liquid crystalline phases at 
all. 

In this paper we take a slightly different perspective. 
We explicitly calculate the orientation-dependent second 
virial coefficient of rather short charged rods numerically, 
for colloidal charges and screening lengths that vary over 
many decades. Such calculations, in which we use ex- 
pansions in spherical harmonics, do not require only the 
asymptotic far-field exp ressions of the multipoles (such as 
considered in Refs. [24l [25l| ) . but in fact their full distance 
dependence. From the resulting second virial coefficient, 



we determine an effective hard-core length and diame- 
ter. Subsequently, we use these — in combination with the 
published hard-core phase diagram [l8[ — to determine 
the expected phase sequence upon increasing the concen- 
tration. This scheme is too crude to distinguish subtleties 
such as whether or not there is a nematic-nematic coex- 
istence regime or to what extent the isotropic-nematic 
phase gap is affected. However, it is supposed to indi- 
cate reliably whether liquid crystalline (N and Sm) or 
plastic crystal (P) phases are to be expected in between 
the isotropic (I) and crystalline (X) phase. We focus 
on the case where the rod length is of the order of the 
screening length or smaller, in contrast to most of the 
previous theoretical work. This is the regime where the 
crossover from N and Sm to P is expected to occur. In 
the limit where the rod length is small and the hard- 
core interactions are important, we give a simplified the- 
oretical description that turns out to be in remarkable 
agreement with the numerical results. As our numerical 
approach relies on an expansion in spherical harmonics 
of the effective pair interaction between two rods, it leads 
to explicit but involved expressions. We present some of 
the mathematical technicalities of the derivation of these 
expressions in the appendix. 

II. MODEL 

We consider a system of identical charged colloidal 
rods suspended in an electrolyte solvent of dielectric con- 
stant e, Debye screening length and Bjerrum length 
Ib = e 2 / (ATreksT), at temperature T. Here e is the ele- 
mentary charge, and fee is the Boltzmann constant. The 
rods are assumed to have the shape of a spherocylin- 
der consisting of a cylinder of length L and diameter D 
capped by two hemispheres also of diameter D. The rods 
have a fixed charge, which we treat here as an (effective) 
line-charge density eX distributed homogeneously on the 
axis of the cylinder. We are interested in the effective 
pair potential V(r; u, &') between two rods with orienta- 
tions u> and Cj' at a center-to-center vector r, thermally 
averaged over the degrees of freedom of the electrolyte 
solvent (characterized by k _1 and Ib)- In the spirit of 
Derjaguin, Landau, Verwey and Overbeek (DLVO), we 
assume that the effective pair potential consists of steric 
hard-core repulsions and electrostatic screened- Coulomb 
interactions between segments of the line charge of the 
two rods. We ignore short-ranged Van der Waals attrac- 
tions (i.e. we assume the particle and the solvent to be 
index-matched or that the dispersion forces are cancelled 
by steric or charge stabilization). Within these approxi- 
mations the effective pair potential can be written as 

{oo for overlapping rods, 
(1) 
/3V C (r; &,&') otherwise, 

where = /cbT 1 , the overlap refers to the hard-core 
repulsions, and the electrostatic interaction potential is 
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given by 

0V; (r; <&,£>') 



dl / dZ 



, exp[— n\r + l'u>' — lu\] 
\r + l'Co' -lCo\ ' 

(2) 

The integration variables I and V play the role of coordi- 
nates running along the cylinder axis of each of the two 
rods, from one end of the cylinder to the other end. In the 
long- rod limit, L/D ^> 1 and kL ^S> 1, one can replace 
the integration domains in equation @ by the full real 
axis, together with the constraint that the cylinder axes 
are in "cross configuration" (i.e. the axes intersect when 
projected onto the plane parallel to both axes). Other- 
wise, the potential vanishes. One then easily shows that 
V e only depends on the shortest distance and the relative 
angle between the two rods [1, HH H3] • Here we focus on 
shorter rods, for which this simplification does not ap- 
ply. In the appendices we derive systematic series expan- 
sions in spherical harmonics to describe the angular and 
position dependence of V c explicitly, focussing on rods 
that are rather short compared to the Debye screening 
length (which sets the range of the electrostatic repul- 
sions). More specifically, the expansion of the angular 
dependence is truncated and we consider each term as 
an expansion in kL up to fourth order (see appendix). 
We compare the result with the large- kL limit. 

The present model can be characterized by a few di- 
mensionless combinations. In the limit of uncharged rods 
(A = 0), the aspect ratio L/D of the hard-core dimen- 
sions is of primary importance. However, for the charged 
rods of present interest, the ratio kL (of the hard-core 
length to the Debye screening length of the solvent) gives 
more information on the interaction anisotropy. The ra- 
tio kD is relevant as a measure of ionic strength. Dimen- 
sional inspection of the expression in equation (J2]) shows 
that the strength of the electrostatic interactions is de- 
termined by the dimensionless (square of the) line charge 
density 



(3) 



These dimensionless combinations can span quite a range 
of numerical values in experimental systems. For in- 
stance, for fd virus suspended in water one finds [|[ 
L/D > 100, kD ~ 0.1-1, and q = 70-700, and recently 
synthesized silica dumbbells in oily solvents [4l| are best 
characterized by L/D ~ 1, kD ~ 1, and q ~ 10 2 . Short 
(double stranded) DNA chains have kD ~ 0.1-1 and 
q ~ 0.1-10, while their length can be varied by the num- 
ber of base pairs included in the sequence. These chains 
can be characterized as rigid rods up to the persistence 
length corresponding to L/D ~ 50. Moreover, present- 
day synthesis techniques allow for the tuning of surface 
charge, in principle at least, from essentially vanishing 
to extremely high. This is achieved for example by using 



different coatings with varying degrees of ion-dissociation 
of the surface groups. It is therefore of interest to inves- 
tigate the thermodynamics of the present model over a 
wide range of parameters. 



III. THERMODYNAMICS AND EFFECTIVE 
DIMENSIONS 

With the pair potential specified by equations ((TJ and 
Pj), and with an explicit scheme to evaluate it as ex- 
plained in the appendix, we can study the macroscopic 
properties of suspensions of these charged rods. In princi- 
ple, we do this as a function of concentration, for various 
q, kD, and L/D. Here we circumvent the complexity 
of the full statistical-mechanical calculation of free en- 
ergies and phase diagrams of the system at hand. We 
do this by mapping the second virial coefficient of the 
charged spherocylinders of interest onto that of hard sphe- 
rocylinders with an effective cylinder length L e g and an 
effective diameter £> e ff that we will calculate below. We 
then presume that the phase diagram of the system of 
charged rods follows from that of the effective hard-rod 
system, which we take from published computer simula- 
tion data [HI . It is well-known from these and follow-up 
simulations of hard rods, as well as density functional 
theory [26|, [27j , that this system exhibits a sequence of 
phase transitions upon increasing the concentration that 
strongly depends on the aspect ratio L/D: sufficiently 
elongated hard rods with L/D > 3.7 have a phase se- 
quence isotropic-nematic-smectic-crystal (I-N-Sm-X) , 
sufficiently short hard rods < L/D < 0.35 show a se- 
quence I-P-X with P a plastic crystal, and in between 
there are two more regimes in which the N and P phase, 
respectively, do no longer appear in the phase sequence. 
Below we determine how the analogous crossovers be- 
tween these regimes of the effective system, as determined 
by L e g/D e f{, depend on the parameters q, kD, and L/D. 

A key ingredient of our calculation is the effective ex- 
cluded volume E(lo,lo') of two charged rods with orien- 
tations Co and Co' defined as 



E(ui,ui') = J dr (l-exp[-pV(T;w,w')] 



(4) 



where the pair potential between the rods is given in 
Eqs. ([T]) and Note that E(Co,Co') is in fact twice 

the corresponding second virial coefficient, and that the 
nomenclature "effective excluded volume" stems from the 
fact that it reduces to the actual excluded volume of the 
pair in the case of purely hard-core interactions. On the 
basis of symmetry arguments one easily checks that the 
angular dependence of E(u), Co') is in fact only through the 
angle 7 = arccos(u) • Co') between the cylinder axes of the 
two rods. In Fig. [1] we show this 7-dependence of E for 
rods characterized by kL = 1 and kD = 0.01 (so L/D = 
100 and weak screening), for several charge parameters 
q ranging from q = (uncharged) to q = 0.01 (fairly 
charged) . The results of Fig. [1] stem from a combination 
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tq 0.08 
0.06 




Figure 1: The effective excluded volume k 3 E, as a function 
of the angle 7 between the two rod orientations, for different 
values of the charge parameter q. We used the parameter 
values kL = 1 and kD = 0.01 (such that L/D = 100). 



of numerical and analytic procedures explained in detail 
in the appendix. These involve a five-fold integration: 
over the contour of the rods I and I' in Eq. ([2]), and the 
center-to-center separation vector r in Eq. (HJ) . 

The key observations of Fig. [TJ which is typical for 
many system parameters, are that for increasing q the ef- 
fective excluded volume becomes (i) less anisotropic, and 
(ii) larger in magnitude. Moreover, for all q the effective 
excluded volume is larger for perpendicular orientations 
than for parallel ones. Qualitatively, and in fact quanti- 
tatively for many parameters, this behaviour is identical 
to that of hard spherocylinders of effective length L e g- and 
diameter D c g, for which the excluded volume is given by 



3] 



4l7T 



2irL cS D 2 cS + 2L 2 cS D eS sin 7. (5) 



In principle one can fit the functional form of Eq. ([S]) to 
the numerical results such as those of Fig. [1] to deter- 
mine the effective hard-core dimensions L e ff and D c ff for 
given charged-rod parameters. However, instead of fit- 
ting the full angular dependence numerically, it is more 
convenient to match the isotropically-averaged effective 
excluded volume and the parallel one, given by 



E\c a — 



1 



(4*) 2 J dd> J d "' E ^ 

1 r 

- I d 7 sin 7^(7), (6) 

v U = E(u,u) =£(7 = 0), (7) 

to the values for spherocylinders with effective hard-core 
dimensions L e ff and D c ff 

2irL cB Dl s 

47T 



Viso — —D eS 



2ttL c bD cS , 



(8) 
(9) 



respectively. This procedure yields the effective hard- 
core dimensions 



3£ 



II 

47T 



+ 3A - V3A(2 + 3A)) 



^ = 2A + ^3A(2 + 3A), 



(10) 

(11) 



where we used, for notational convenience, the dimen- 
sionless anisotropy parameter A denned as 



A = 



E 



En 



(12) 



It turns out that inserting L g and D G ff as obtained from 
Eqs. ([TO]). (fTTj) . and (HU) into Eq. © gives an angular 
dependence that is in very good agreement with the nu- 
merically obtained effective excluded volume of charged 
rods. 

It is also interesting to compare our numerical results 
with analytic expressions that are valid in the limit where 
L/D 3> 1 and nL » 1, as obtained by Stroobants et al. 
}3l| . In this needle-limit the effective excluded volume is 
given by 



£^(7) =2L 2 K - 1 sin7 



7e + In 2nq — In sin 7 
2nq exp[— kD] 



r 



sin 7 



, (13) 



where 7e ~ 0.577 is the Euler-Mascheroni constant and 
where the incomplete gamma function (or exponential 
integral) is defined by 



T(a, x) 



/>oo 

/ dyy"' 1 exp[-y]. 

J X 



(14) 



From this expression — using the Onsager limit Vi so ,oo = 
(7r/2)L 2 D c ff,oo for the isotropically averaged excluded 
volume — the effective diameter can be calculated 

K-Doff,oo = 7e + In 27rg + In 2 - - 

'd7sm 2 7 rfo, 27r9eXp[ - KjP] V (15) 

The effective length is taken equal to the rod length 



IV. NUMERICAL RESULTS 

Calculations such as those of Fig. [T] are reasonably ac- 
curate for values of kL roughly up to 2. For higher values 
the applied approximations become poor, such that for 
kL > 3 the calculations become even qualitatively un- 
reliable for many of our parameters. For this reason we 
restrict most of our attention to the regime where kL < 2. 

In order to assess the accuracy of our calculations, 
we compare some of the results of our calculations with 
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Figure 2: The effective diameter K_D e ff as a function of the rod diameter kD for (a) reL = 1 and different values for the charge 
parameter q; (b) q = 10 and different values for the rod length kL. The rod dimensions are scaled by the screening length «~ . 
The thin solid line is a guide to the eye, representing the hard-core limit D B g = D. The small solid circles give the values for 
kD — from the numerical calculations. The larger open circles are obtained by the approximation given in Eq. (|17[) . 



those obtained from more extensive numerical integration 
schemes. One is given by the same spatial integration 
scheme as before, but with the (effective) line-charge den- 
sity replaced by a discrete charge distribution. The rod 
charge is represented by an odd number of charge units 
(AI = 2N + 1) distributed evenly on the cylinder axis, 
where one unit is always located on the center of the axis, 
and two units are always located on the two end points 
of the axis. The latter are of magnitude eXL/(4N), while 
all others are of magnitude eXL/ (2N). This ensures that 
the total charge is eXL and the continuum limit N — > oo 
yields the correct homogeneous line charge. The other 
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Figure 3: The effective excluded volume k 3 E, as a function 
of the angle 7 between the rod orientations, calculated us- 
ing different numerical schemes (see text), involving M dis- 
crete charges, Monte-Carlo (MC) integration, and the present 
analytic approach. We used the parameter values kL = 1, 
kD = 0.25 (such that L/D = 4) and q = 1. 



scheme uses the same discrete charge density as described 
above, but uses a Monte-Carlo (MC) scheme to perform 
the integration. This scheme is denoted by the plusses in 
Fig. [3J The agreement between the results obtained from 
the different schemes, as shown in Fig. [3J is excellent for 
M > 13, particularly when considering that the shape 
of the effective excluded volume differs significantly from 
the hard-core case for these parameters. Therefore, we 
conclude that our calculation correctly predicts the angu- 
lar dependence of the effective excluded volume of short 
charged rods. 

In the previous section, we have shown that the an- 
gular dependence of the effective excluded volume can 
be used to calculate the effective rod dimensions L e g 
and -D e ff — from the values of E\\ and E± — by applying 
Eqs. COB, CP, and |(T2J). Fig.^a) shows the numerically 
calculated effective diameter as a function of the real di- 
ameter for kL = 1 and a range of charge parameters q. 
Fig. HJb) shows the same function, but then for q = 10 
and a range of rod lengths kL. Note that all (effective) 
rod dimensions are expressed in units of the screening 
length. In Fig. [2jb) the needle limit kL ^> 1, given by 
Eq. (fT5")) . is plotted for comparison. Both graphs clearly 
reveal two regimes 



D e for D<^D C 



D for D > D e 



(16) 



These can be identified as an electrostatic regime at 
small kD (weak screening) and a hard-core regime at 
high enough kD (strong screening). In the hard-core 
regime, the effective diameter equals the hard-core diam- 
eter, while in the (weakly screened) electrostatic regime 
the effective diameter saturates to a plateau value D e . 
This electrostatic effective diameter depends on q and 
kL, and increases with increasing q and kL. Also, it is 
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Figure 4: The effective length L g / L as a function of the rod diameter kD for (a) kL — 1 and different values for the charge 



parameter q; (b) q = 10 and different values for the rod length kL. 
the effective length is scaled by the rod length L. 



The rod diameter is scaled by the screening length k and 



(much) larger than the hard-core diameter due to the 
(strong) rod-rod repulsions. Values of the electrostatic 
effective diameter are included in FigfSl where the small 
solid circles represent values obtained from numerical cal- 
culations for kD = 0. The larger open circles represent 
the following simple approximation for D e . 

In the short-rod limit, we can treat the double layer 
around the rod as spherically symmetric, with an effec- 
tive point charge eXL in the center, such that also the pair 
potential is spherically symmetric. This gives A = 0, and 
hence from Eq. (fTX)]) . for large enough q (or small enough 
kD), we obtain the electrostatic effective diameter from 
the simple expression 



dx x 1 — exp 



-qn L 



2 r 2 ex P[ 



This approximation is given in Fig. [5] by the larger open 
circles. Both graphs show good agreement for kL < 1 
and all values for q. Fig.[5Jb) also shows that the regime 
kL < 2 — which is reliably accessible with our truncated 
numerical scheme — evolves smoothly to the needle-limit 
kL ^> 1 of Stroobants et al. [3l|. The curve for kL = 3 
shows some signatures of the numerical instabilities we 
encounter for larger kL. 

In a similar fashion wc can also study the effective 
length L c ff of the rods. Fig. SJa) shows results of numer- 
ical calculations of the effective rod length for kL = 1 and 
a range of charge parameters q. Fig.^b) is the result for 
q = 10 and a range of rod lengths nL. The rod dimen- 
sions are expressed in units of the Debye length, whereas 
L e s is expressed in units of the hard-core length. We dis- 
tinguish again two asymptotic regimes, the strong screen- 
ing (hard core) regime kD 3> 1 where L c ff — L, and the 
weak-screening (electrostatic) regime kD <C 1 where L e ff 
reaches a plateau value that depends on q and kL. Note 
also that L c g < L which is perhaps unexpected at first 



sight. Naively, one could expect the effective length to 
increase with increasing effective excluded volume. How- 
ever, as Sato and Teramoto [34[ pointed out, the effective 
length decreases with increasing rod charge density be- 
cause of end effects. Thus, the increase of the effective 
excluded volume — due to the increase of the rod charge 
density — is purely caused by the increase of the effective 
diameter. Moreover, this increase balances the decreas- 
ing in effective length such that the total effective par- 
ticle length L c g + D e g does increase with increasing rod 
charge density. Inspection of Fig. IDJa) also reveals nu- 
merical (convergence) problems for q > 100 at kD > 1, 
where nL c ff sharply drops and rises before reaching the 
hard-core limit L e s — L. This is in fact only a mi- 
nor problem in practice, as it only occurs in the regime 
where L e e/D a g < 0.1. There, the anisotropic contri- 
bution to the effective excluded volume is much smaller 
than the isotropic part. Upon approach of the needle- 
limit kL 1, see Fig. |U[b), we find that L g approaches 
L for all values of kD, as expected. 



V. PHASE BEHAVIOUR 

We have determined the effective length and diam- 
eter of charged rods, by mapping their orientation- 
dependent second virial coefficient onto that of effec- 
tive hard rods. Subsequently, we also study the effec- 
tive length-to-diameter ratio L e ff/D c g. In Fig. [5] we 
show this effective aspect ratio as a function of the rod 
charge for kL = 1 and a range of rod diameters kD. All 
curves with kD > essentially decrease from their max- 
imum value — the hard-core aspect ratio L/D — towards 
the curve given by kD = 0. This indicates that the ef- 
fective dimensions of charged rods become independent 
of the hard-core diameter for large charge parameters, 
where we enter the electrostatic regime. Also, since the 
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Figure 5: The effective aspect ratio L e g/D e g as a function of 
the charge parameter q for kL = 1 and different values for 
the rod diameter kD. Different — possibly coexisting — phases 
are associated with a certain range of (effective) aspect ratios. 
See the text for an explanation of the abbreviations and the 
boundary values. 



effective aspect ratio for kD = is a decreasing func- 
tion for large q, we see that the charged rods essentially 
behave like charged spheres upon increasing the charge 
above a certain value. 

Moreover, Fig. [5] reveals a local maximum for very 
small kD, in the regime where q ~ 1. This effect can 
be understood by considering the electrostatic regime 
for small charge parameters q. Eq. (fTTj) shows that 
the effective aspect ratio is governed by the dimcn- 
sionless anisotropy parameter A, which is defined in 
Eq. (fT2"| . In the electrostatic regime, this anisotropy 
can be shown — up to first order — to be proportional to 
q. The reason for this is that the linear approxima- 
tion of the effective excluded volume is orientation in- 
dependent ■ Therefore, the difference between the 
isotropically-averaged and parallel values is of second or- 
der in q, whereas the parallel value itself is of first order. 
The effective aspect ratio is of order \fK, and thus in- 
creases as y/q. Conversely, for q > 1 the effective length 
is more or less constant, and the effective aspect ratio 
decreases again due to the increase of the effective diam- 
eter. 

The horizontal dotted lines in Fig. [5] indicate the 
crossover values (0.35, 3.5, and 3.7) for regimes with dif- 
ferent phase sequences. The values for these aspect ratios 
are taken from simulation results of hard-spherocylinder 
systems by Bolhuis and Frenkel 18]. These simulations 
consist of explicit free-energy calculations of coexisting 
phases, where the most dilute phase is always given by 
an isotropic fluid (I), and the most dense phase by a fully 
ordered crystal (X). Depending on the aspect ratio, dif- 
ferent phases were found in between these two phases. 
For aspect ratios exceeding ~ 3.7 the phase sequence I- 



N-Sm-X was found upon increasing the density. Here, 
the N and Sm denote the nematic and smectic-A liquid 
crystalline phases, respectively. Somewhat shorter rods, 
with an aspect ratio in the narrow regime between ~ 3.5 
and ~ 3.7, can still form a smectic-A but no longer a ne- 
matic phase, and hence have a phase sequence I-Sm-X. 
Even shorter hard rods, with an aspect ratio in between 
~ 0.35 and ~ 3.5 cannot form a thermodynamically sta- 
ble smectic-A phase, and thus crystallize directly into a 
fully ordered crystal from the isotropic fluid, yielding a 
phase sequence I-X. Very short hard rods, with an as- 
pect ratio smaller than ~ 0.35, exhibit a plastic (P) crys- 
tal phase, such that the phase sequence is I-P-X. The 
plastic crystal phase is characterized by orientational dis- 
order, but has translational order as in a crystal phase 
[lil I33 ] . This regime arises naturally in the case that kL 
is small. Then, such a crystal forms because of the essen- 
tially isotropic long-range repulsive interactions, but the 
competition with entropic effects prevents the rods from 
aligning. 

We use the mapping of the charged-rod system onto 
the effective hard-rod system to give an indication of the 
phase sequence of systems of charged rods as a function 
of the parameters kL, kD (or L/D), and q. For instance, 
from the curve for kD = 1 in Fig. we see that the ef- 
fective aspect ratio never exceeds unity for any q. This 
excludes the possibility of a nematic or smectic-A liquid 
crystal phase. The curve starts off at its maximum (in 
the limit where q — ► 0), where the effective aspect ra- 
tio equals the hard-core aspect ratio L/D = kL = 1. 
It crosses the value L c g/D c g — 0.35 at q « 2.35, such 
that a sufficiently large rod charge density allows for a 
plastic crystal phase. Similarly, for kD = 0.1 (which cor- 
responds to L/D = 10), we find all four phase sequences 
upon increasing q. 

By determining the intersections of the effective aspect 
ratio with the crossover values of the hard-rod system, 
we construct "phase diagrams" indicating the different 
regimes. In Fig. [5] we present two examples of such dia- 
grams in the plane spanned by q and kD. In Fig. [Ha), 
we fix kL = 1, such that the horizontal axis could read 
D/L as well. In Fig. E^b) we fix L/D = 20, such that 
the change in kD physically corresponds to a change in 
salt concentration (while keeping the particle dimensions 
fixed). The symbols denote the crossover values for the 
effective aspect ratio as determined from our numerical 
data (such as presented in Fig. [5]). The lines are based 
on an approximate theoretical model to be discussed in 
section IVII 

Both diagrams in Fig. [5] show that rods with suffi- 
ciently high surface charge density always show the I— 
P-X sequence. This is due to the essentially spherical 
nature of the effective shape of highly charged rods. The 
limit of uncharged rods is determined by the hard-core se- 
quence that corresponds to L/D. The I-N-Sm-X regime 
at fixed kL in Fig. is completely bounded. First, 
by a hard-core regime when kD > 0.27, where the liq- 
uid crystal phases cannot exist even for q = because 
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Figure 6: Boundary lines for given values for the effective aspect ratio L e g/D e g. See the text for an explanation of the 
abbreviations of the different regime labels. The points are results of the numerical calculations, and the lines are given by a 
simplified theory. We fix (a) nL = 1, and (b) L/D — 20, respectively. 



L/D < 3.7. Second, by an electrostatic regime in the 
weak-screening limit of small kD, where the rods effec- 
tively behave as spheres since D c g 3> L c g. Conversely, 
the trends displayed for fixed L/D in Fig.JSJb) are mono- 
tonic, with an I-N-Sm-X regime that extends to higher 
q with increasing kD. 



effective charge on the rods 

- 4-7T - 

m = y 5( 7 ) 3 

r°° , / 

+ Air drr 1 — exp 

Jd( 7 ) \ 



-qn 



2 r 2 CX Ph CT ] 



KT 



(20) 



VI. A SIMPLER MODEL 

For small values of the effective surface-charge density, 
we found that the electrostatic contribution to the effec- 
tive excluded volume is essentially isotropic in nature. 
This means that the anisotropic effects are primarily due 
to the hard-core anisotropy (as apparent from Fig. [I}, 
such that 



L D. 



(18) 



On this basis, we propose here a simple model, which 
turns out to describe our numerical findings with re- 
markable accuracy. This model introduces a "spherical 
approximation" of the electrostatic contribution to the 
effective excluded volume, which involves the orientation- 
dependent diameter D(j). The volume of a sphere of this 
diameter is equal to the hard-core excluded volume of a 
pair of rods 



47T - . . , 47T , 

—DM 3 = —D 3 

3 K1J 3 



2ttLD 2 + 2L 2 D sin 7. 



(19) 



We approximate the effective excluded volume by the 
value for a charged sphere of diameter D (7) and an ef- 
fective surface charge that equals the total amount of 



Note that the only orientation dependence of the elec- 
trostatic contribution to this effective excluded volume 
(i.e. the second term) comes from the integral boundary 
D(fy). To calculate the effective dimensions, we only need 
the parallel and isotropically averaged values of the ef- 
fective excluded volume. In the parallel case (7 = 0) this 
value is readily calculated 



47T _, 



- 47T 



where 



dr r 2 



47T 



1 — exp 



4-7T 



9 9 exp[— nr] 
-qn z L z —^ 



(21) 



—Di = -^D 



2ttLD z . 



(22) 



The isotropically-averaged value can be calculated nu- 
merically by using expression (|20p . However, we approx- 
imate it by the value for a charged sphere of diameter Z?i so 
(using the same total effective charge), which is taken 
from the isotropic average of the hard-core excluded vol- 
ume 



— D 3 

3 lso 



47T 



-D 3 + 2txLD 2 + —L 2 D. (23) 
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Figure 7: (a) The effective length nL c g as a function of the rod length kL for nD = 0.1 and different values for the charge 
parameter q. The thin solid line represents the needle or hard-core limit, where L e g ~ L. (b) The effective diameter nD a g as a 
function of the rod charge parameter q for kD = 0.1 and different values for the rod length kL. The (effective) rod dimensions 
are scaled by the screening length kT 1 . 



This approximation yields the simple expression 



- 4-7T - 3 

° = ~3~ 80 
+ Air dr r [ 1 — exp 



9 9 exp[— kt] 
-qtCL 1 — — - 



(24) 



With our explicit expressions (|2T|) and (|24j) . we eval- 
uate the effective dimensions from Eqs. (fT0|) and (fTT|) 
as before. The resulting crossover values of the hard- 
rod system are shown by the curves in Fig. [6l and are 
in very good agreement with the numerical calculations 
(denoted by the symbols). The key to this remarkable 
accuracy lies in the fact that the anisotropic electro- 
static contributions are relatively unimportant, because 
the rod length is small with respect to the screening 
length (i.e. kL < 2). Thus, our simple model accounts 
for the hard-core anisotropy correctly, as well as for the 
isotropic electrostatic contribution. 

In a sense, this theoretical description can be viewed 
as a kind of perturbation theory, where we expand the 
pair potential as a function of kL. The hard-core re- 
pulsion represents the zero-th order. The lowest-order 
contribution to V c is quadratic in kL and independent 
of rod orientations. Also, it happens to correspond to 
the interaction potential of two point charges eXL. If we 
plug this approximation of V(r; u>, a)') into the expression 
of the effective excluded volume (given by Eq. ([4|), we 
obtain an expression where the integral boundary D is 
still a function of both the angle between the rod orienta- 
tions and the direction of the center-to-center separation 
vector r. In fact, it is given by the distance where the 
rods touch, given a certain orientational configuration. 
By setting this overlap diameter to a value that is inde- 
pendent of the orientation of r, but still respects the total 
hard-core excluded volume, we effectively neglect its de- 
pendence on kL. This choice is justified by the fact that 



(for small kL) the size of the double layer around the par- 
ticles is larger than the variations in the overlap diameter 
D. That is why our simple theoretical description can be 
interpreted as a perturbation theory of a hard-rod refer- 
ence system with an (almost) isotropic electrostatic con- 
tribution. Unfortunately, it completely fails to describe 
the anisotropic effects in the electrostatic regime. In this 
regime the anisotropic details of the electrostatic contri- 
butions do become important compared to the hard-core 
contributions. 



VII. DISCUSSION AND CONCLUSION 

The numerical results presented in this paper give ac- 
cess to a part of the parameter space where there is a 
large difference between the effective length and the real 
length. In this regime, one cannot hope that the theory 
of Stroobants et al. [3l| gives any accurate results, as this 
is based on the needle limit where L^a ~ L. The pertur- 
bation theory of Chen and Koch [37| breaks down for 
most of our parameter values. This is because it is based 
on small charges, and thus fails to describe the effect of 
large rod surface-charge densities. Also, this theory is 
not accurate for large differences between the effective 
and hard-core diameter. 

In Fig. [7(a) , we show results of numerical calculations 
of the effective rod length as a function of the hard-core 
length, for kD = 0.1. Note that again the effective length 
is always smaller than (or equal to) the hard-core length. 
Also, in accordance with the results from Fig. HJ there is 
a hard-core regime for small values of the charge parame- 
ter q, as well as for small values of the rod length kL, for 
which the total amount of effective rod charge is small. 
On the other hand, there is an electrostatic regime. In 
Fig. [U this was shown to be the case for decreasing val- 
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ues of kD, where the plateau value (i.e. the electrostatic 
length) depends on q and kL. However, from Fig.[7{a), it 
can be seen that this electrostatic length depends mostly 
on the rod length kL, and not really on the charge pa- 
rameter q, as long as either q or kL is large enough. Fur- 
thermore, the effective length is "wedged" in between the 
electrostatic length and the hard-core length, where the 
electrostatic length approaches the hard-core length in 
the needle limit (kL 3> 1). Unfortunately, there is no an- 
alytic theory yet that describes our numerical results for 
this electrostatic length as a function of kL. Therefore, 
it would be worthwhile to gain new insight in the effect 
of electrostatics on the effective rod length for interme- 
diate kL — neglecting hard-core interactions — in the case 
of large rod charges. Additionally, Fig. Efb) shows re- 
sults of numerical calculations of the effective diameter 
as a function of the charge parameter q. For q£l, there 
is a smooth transition to the theoretical needle limit of 
Ref. [3l|, where kL ^> 1. Conversely, this is not the case 
for q < 1, due to the fact that the approximations lead- 
ing to Eq. (fT5|) do not give the correct effective excluded 
volume for small values of q and (nearly) parallel rods. 
More investigations need to be made into this regime. 

In conclusion, we have numerically studied the second 
virial coefficient of short charged rods dispersed in an 
electrolyte, presuming pairwise screened- Coulomb inter- 
actions between the line-charge segments of the rods. 
The control parameters of interest are the hard-core 
length L and diameter D, the Debye screening length of 
the medium and the charge parameter q. The main 
resulting quantities are the effective diameter D e g and 
length L e g of the rods. By a mapping onto an effective 
hard-core system — for which the sequence of phases be- 
tween the dilute isotropic phase and the dense crystalline 
phase is known for all aspect ratios — we predict the rela- 
tions between control parameters and the expected phase 
sequence explicitly. We have also constructed a simplified 
model, based on the diameter D(^) of Eq. (fl§j) , which re- 
produces the numerical results accurately at the expense 
of much less computational effort. This model is partic- 
ularly successful in the regime of large effective aspect 
ratios (L e s/D c g > 1) and small ratios of the rod length 
to the screening length [kL < 1). 

An important result of this work is that highly charged 
short rods at low salt concentrations (i.e. at strong 
Coulomb couplings) have a strong tendency to form plas- 
tic crystals upon compression. The plasticity stems from 
the large effective diameter, which make the rods behave 
essentially as inflated repulsive spheres with only small 
nonspherical interactions that are too weak to cause ori- 
entational ordering in the crystalline phase. This finding 
could be important in the study of silica or gold nanorods, 
that have reasonably large hard-core aspect ratio (like 
L/D ~ 5). Here, liquid crystalline phases could be ex- 
pected, but only if the charge on the rods is small enough. 
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Appendix A: THE PAIR INTERACTION OF TWO 
CHARGED RODS 

The pair interaction of two charged rods is given by 
Eq. |(3J), where we assume that the electrostatic inter- 
action is determined by integrating over pairs of effec- 
tive line-charge elements interacting with the screened 
Coulomb potential. The distance between these pairs is 
given by a superposition of the relative position of the 
rods and the combination of the position of the line el- 
ements along both rods. Since the integral in Eq. ([2]) 
cannot be calculated analytically, we try to simplify the 
calculation. By expanding the integrand in spherical har- 
monics, we obtain terms that factorize into two functions 
of the respective positions 

eXP , [ ~ |r 7 8|] - V>2 + l)h(r)Pi(r ■ for r > s, 

|r-s| ^ 

oo +i 

1=0 m=-l 

(Al) 

where ii and ki are the modified spherical Bessel func- 
tions of the first and second kind, respectively. These 
functions are given by 

= ^I l+h {x), (A2) 

k l (x) = x P^K l+h (x), (A3) 

where /„ and K v are the modified (cylindrical) Bessel 
functions of the first and second kind, respectively. The 
Legendre polynomials Pi are expanded into spherical har- 
monics Yi. m using the famous addition theorem. We use 
the notation where r — |r|, and f = r/r. Finally, the 
asterisk "*" denotes complex conjugation. The unit vec- 
tor as given in the arguments of each of the spherical 
harmonic functions should be interpreted as the two an- 
gles in spherical coordinates with respect to an arbitrarily 
chosen reference frame. Since the Legendre polynomials 
of the dot product of the two orientations is independent 
of this choice, so is the sum over to of the product of the 
two spherical harmonics. 

We note that one could consider rewriting the expres- 
sion of the pair potential in rotational invariants (as used 
in Ref. [4fJ). These are functions of three orientations, 
including a sum over to of a product of three spherical 
harmonic functions multiplied by Clebsch-Gordon coeffi- 
cients. They form a complete set of orthogonal functions 
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dependent only on the relative orientations of f , w, and 
u>' with respect to each other . However, it turns out 
that in our case these are not really helpful. Alterna- 
tively, one could consider a resummation of the expan- 
sion in spherical harmonics, such that each term has a 
faster asymptotic decay than the previous term. This is 
not the case here, since each Bessel function fc; has the 
same asymptotic decay as ko [24j . 



Appendix B: DOMAINS OF INTEGRATION 




Figure 8: Illustration of the domain of integration of the su- 
perposition of the positions of two line elements. The dashed 
circle of radius r divides the parallelogram into two domains. 

The integration over line elements of both rods in 
Eq. @ is in fact an integration of the vector ICj — I'd' 
over a parallelogram-shaped area in the plane tangent to 
both rod orientations. This area is illustrated in Fig. [8] 
There is a straightforward choice for the reference frame 
and a substitution of variables 



7 7 N 
cos-, sin-, 0^ 



' cos — , — sin — , 



2 2 
ICj — I Q) — {p cos p, p sin ip, 0) , 



(Bl) 

(B2) 
(B3) 



where 7 is the angle between the two rod orientations. 
The polar coordinates p and ip describe the same plane 
as I and I'. The parallelogram can be cut up into four 
equivalent pieces, keeping only the terms in the expan- 
sion (|A1[) where / and to are both even. The integral 
boundaries of the first quadrant (0 < <p < tt/2) satisfy 



< p < 



L sin 7 



2 sin (<p+i) 



(B4) 



It is important to note that the functional form of the 
integrant can vary as a function of p, because ki and ii 
switch roles when r < p. We shall split the result of our 
expansion into each order in I and to, to be examined 



separately. We write 



x Ai, m (r;>y)Pi, m (cos9) cos(to.(/>), 

(B5) 

where P\^ m are the associated Legendre functions. We 
have used that for I and to both even 



+ Yl,_ m (M) 



<2Z + 1 (/ 



4-7T (7+ to)! 



Pi,m (cos 6) cos(m</>) , 



(B6) 



and 



2 ( Y[ * m ( 



IT 



Y, 



l,-r 



IT 



21 + 1 yj (l + m)\(l - m)\ 



4tt 2 l ( l+ ™ ) ! ( '~ 2 m )! 
The integral Ai : m(r;j) in Eq. (|F35|) is given by 
4 



cos(mip). (B7) 



-4/,m(V;7) 



sm 7 J 



dip cos(mip)Bi (r ; ip, 7) , (B8) 



where 



L Bin y 

Bi(r;<p,j) = ki(nr) / 2 " n ( v+ 2) dppii(np) 
Jo 



(B9) 



for r > — L / 1i17 ^n , and 

2sin( V +i) ' 



Bi(r;ip,j) = ki(nr) dppii(np) 
Jo 



+ ii(nr) 



dppk(Kp) (BIO) 



for r < L f n ' 1 , l , . 

Let us have another look at Fig. [5] The dashed circle 
indicates the value for which the variables r and s in 
Eq. (|A1|) switch (in this case s is replaced by p) . Consider 
the first quadrant (i.e. the upper right-hand corner). Let 
us also assume 7 < tt/2. In the end, we will calculate 
the effective excluded volume for < 7 < it, but this 
expression is symmetric in 7 «-* tt — 7 (due to up-down 
symmetry) so we need only the first half of this interval. 
We describe the integral boundary for p as a function of ip 
just as we describe the boundary of the parallelogram by 
p as a function of ip. However, the integrand in B changes 
when the boundary of the parallelogram intersects with 
the circle of radius r. Therefore — depending on the value 
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of r — we have one to three domains for B as a function 
of <p 

¥>e[0,f] forr<^, 
tpe [0,a(r)],[a(r),/3(r)],[/?(r),|] 

for <r <Lsin^, 

93 € [0, a(r)] , [a(r), f ] for i sin ^ < r < L cos ^, 
<^ G [0, §] for r > Lcos^, 



where 



a(r) = arcsin 



L sin 7 \ 7 



2r 7 2' 



(Bll) 



p(r) = tt - arcsin (^^) - | , (B12) 

are the angles for which the circle intersects the bound- 
ary of the parallelogram. In each domain, we calculate 
the integral Ai >m using the corresponding expression for 
the integrant Bf (|B10|) if the circle segment lies in the 
interior of the parallelogram; (|B9|) if it lies outside of the 
parallelogram. 



Appendix C: THE LIMIT FOR PARALLEL RODS 

In principle, calculations of the effective excluded vol- 
ume for parallel rods involves the limit 7 — > of 
Eqs. (IB8|) (|B12j) . To obtain the correct result, one has 
to take care to perform the limit correctly in each ex- 
pression, which is not straightforward. It is much easier 
to re-evaluate the expressions in this limit analytically, 
starting with Eqs. (|B1[) - (|B3[) . We use the same reference 
frame, but a different substitution of variables 



w = (1,0,0), 

a' = (i,o,o), 

ICj - I'Cj' = (±:r,0,0), 



(CI) 
(C2) 
(C3) 



where x = \l — V\. Now the integration is performed over 
relative positions of two points on a single line. Half of 
the combinations is positive (I > I'), the other half is 
negative (/ < I'). The integration boundaries of either 
set is given by 



< x < L. 



(C4) 



The length over which each combination I, I' is realized, 
for a certain value of x, is given by L — x. In accordance 
with the previous expressions, we define the integral A 
for parallel rods as 



Ai, m {r-n = 0) = 2fcj(«r) / dx(L - x)h(kx) (C5) 
Jo 



for r > L, and 

Ai. m (r;^ = 0) = 2ki(nr) / dx (L - x)ii(nx) 
Jo 

+ 2ii{nr) I dx (L — x)ki(nx) (C6) 

J r 

for r < L. Note that the expressions are independent of 
m. 



Appendix D: NOTATIONS, INTEGRALS AND 
TAYLOR SERIES EXPANSIONS 

In order to calculate the integral A, we first need to cal- 
culate the integral B by performing the integration — over 
the radial coordinate p — in Eqs. (|B9p and (|B10p . Intro- 
ducing the notation 



dx x ki(x), 



li(z) = / dxxii(x), 
ICi(z) 

we can rewrite B as 
K 2 Bi(r;ip,^) 

k,L sin 7 



(Dl) 
(D2) 



1V ' V 2sin(vJ+3) 



for r > - L Tl lV 



h(Kr)li(Kr) + ii(nr)K.i(Kr) 

- iiMKt ( " L , Bi V^ ) for r < , L f*2x\ - 

(D3) 

Unfortunately, there is no (easy) way to write the ex- 
pressions in Eqs. (|D1[) and (|E>2|) explicitly for arbitrary 
However, one can give explicit expressions (necessary for 
our calculations) for I = 0,2,4. First, the Bessel func- 
tions 



i (z) = — A^, (D4) 

. , , (z 2 + 3)sinh(z)-3zcosh(2:) 

i2{z) = — z , (D5) 

i 4 (z) = 

(z 4 + 45z 2 + 105) sinh(z) - (10z 3 + 105 z) cosh(z) 



k (z) 



Hz) 

k 4 (z) 



cxp(— z) 

z 

(z 2 + 3z + 3) exp(-z) 



(D6) 
(D7) 



(D8) 



(z 4 + 10z 3 + 45z 2 + 105z + 105) exp(-z) 



(D9) 
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Next, their integrals 

1i{z 



J 4 (z 
IC (z 



cosh(z) — 1, 
zcosh(z) — 3sinh(z) 



(z 3 + 35z) cosh(z) - (10z 2 + 35) sinh(z) 
~ 3 

cxp(-z), 
(z + 3) exp(-z) 

z 

(z 3 + 10z 2 + 35z + 35) cxp(-z) 



(D10) 
(Dll) 

8 

~ 3' 
(D12) 

(D13) 
(D14) 

(D15) 



Unfortunately, we cannot perform the subsequent inte- 
gration — of the angular coordinate ip — in Eq. (|B8|) ana- 
lytically, when we try to calculate A. Therefore, we use 
the series expansions (for even I) 



l 2n (z) = 2 2n Y / 



(2n + k)\z 2n+2k + 2 



k=0 



(2n + 2k + 2) (An + 2k + l)!fc! ' 



(D16) 



f2"nH 2 



1 (~l) fc (2fc)!z 2 "- 2fc+1 
2^™ (2n - 2k + l)(2ra - k)\k\ 



E 



2 2n ^ (2n + 2k + l)(2ra + k)\(2k)\ 



(2n + k)\z 



|_2rt+2fc+2 



fc=0 



(2n + 2fe + 2)(4n + 2fc + l)!fc!' 



Finally, we define the specific combination 
Ci(nr) = h(Kr)li(nr) + ii(nr)K.i(Kr), 



(D17) 



(D18) 



which turns out to be given by a relatively simple expres- 
sion (for even /) 



C2n(nr) 



(n\) 2 A (-l) fc (2n + 2fc)! 
(2n)\ (n + k)\(n - k)\(nr) 2k + x 



such that 
C (z) 
C 2 (z) 
Ci(z) 



(2 n n*) 2 



1 exp(— z) 

z z 

z 2 - 6 2 (z 2 + 3z + 3) exp(-z) 



(D19) 

(D20) 
(D21) 



z 4 - 20z 2 + 280 



8 (z 4 + 10z 3 + 45z 2 + 105z + 105) exp(-z) 

3 z^ ' 

(D22) 



Note that in each expression the first term cancels the 
divergence of the second term in the limit where z — > 0. 
Hence, this limit is given by 



C2n(0) = S n fi. 



(D23) 



This property is also reflected in the series expan- 
sion — useful for calculations for small nr — given by 



Ci n (K,r) = (-1) 



» (2"n!) 2 
(2n)l 



T 2.^ ' 



2 p ^k+2n+3 j p ^ fc-2«+2 j y 2 

(D24) 

Note that the terms for even k < 2n have vanishing co- 
efficients. 

The limit of parallel rods has a different set of expres- 
sions. Therefore, we define an additional notation 



Ji{z) 
£i{z) 



dx- — -ii(x), 

n Z 

o 

dx ki(x). 



(D25) 
(D26) 



In this way, we split each integral in Eq. (|C6[) in two parts 
K 2 Ai tm (r;j = 0) 



' 2kL ki(Kr)Ji(nL) 



for r > L, 



= < 



, + 2kL ii(Kr)(Ci(nr) ~ Ci(kL)) for r < i. 

(D27) 



Evaluation of these integrals result in slightly more com- 
plicated expressions, when compared to the expressions 
for J and K. in Eqs. (|Dl^ - (jDT5)) 



Jo(z) = shi(z) + 



1 cosh(z) 



(D28) 



1 , . 2 z cosh(z) + 3 sinh(z) 
JT 2 ( Z ) = - -shi(z) - - + ^ ^, (D29) 



Ji(z) = -shi(z) + 



3z 



(3z 3 + 70y) cosh(z) - (5z 2 + 70) sinh(z) 



where 



shi(z) — I dx 



sinh(x) 



(D30) 



(D31) 
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(D32) 



is the hyperbolic sine integral. 

£o(,) = r ( o,,)-^t£), 

3 (3z 3 + 5z 2 + 70y + 70)cxp(-z) 

U(z) - -1 (0, z) — A . 

(D34) 

In principle, one now has the exact solutions for A up to 
I = 4. However, we need the expressions in Eqs. (|D28[) - 
(|D34|) to provide a well defined limit for the parallel rods, 
to use in combination with the expressions for arbitrary 
orientations (i.e. the series expansions in Eqs. (|D16|) . 
(ID17[) . and (|D19[) V Therefore, it will be convenient to 
also have these expressions in the form of a series expan- 
sion 



Jl n (z) 

oo 



(2n + k)\z 2n+2k+1 



k=0 



(2n + 2k + l)(2n + 2k + 2) (An + 2k + l)!fc! ' 



(D35) 



7e - ln(z) 



(-1)' 



- y — 

2 2n k _f^, (2n-2fe)(2n-2fc + l)(2n- k)W. 



(2 n nl) 2 1 
(2n)\ ~z 



2/i 



(-l) fc (2fc)!z 2 "- 2fe 



1 OO 

— Y 



+ 2 2 " J2 



< (2n + 2k)(2n + 2k + l)(2n + jfc)!(2fc)! 
(2n + k)\z 2n+2k+1 



k=0 



(2n + 2k + l)(2n + 2k + 2) (An + 2k + l)!fc! ' 



(D36) 



factor (kL) 2 . We give some examples of the calculated 
expressions for I — and m = 0. We make the distinction 
between four domains in r. For r < L s ™ 7 



(D33) K 2 A , Q (r;j) 



dip Cq(kt) - i (nr)ICo 



sin 7 Jo 
27r ( 1 exp(— nr) 



kL sin 7 
2sin(W |) 



sin 7 \ «;r 

sinh(«x) 9 9 
— - k L 



27r sinh(«x) 
sin 7 Kr 



Inltan^tan^ 7 



4 4 

sin 2 7 K 3 i a sin 4 7 



+ \/l + sin 7(2 — sin 7) 



kL 
kL 
12 



24 



2560 



^3 L 3 



y/l + sin 7 (16 - 8 sin 7 + 2 sin 2 7-3 sin 3 7) — 



1 — (7 + 5 cos 2 7) 



3G 



21600 



(El) 



Appendix E: TRUNCATION AND SOME 
EXAMPLES OF EXPRESSIONS 



In principle, the calculation of each of the terms in 
Eq. (|B5j) (i.e. each order of I and m) involves an infinite 
series expansion in kL. We will restrict our calculations 
to / = 0, 2, and 4, and truncate each series expansion. 
Since the integration domain of A is shaped like a par- 
allelogram with sides of length L, we divide out a factor 
L 2 to make both A and B dimensionless (i.e. we calculate 
k 2 A/(kL) 2 and k 2 B/(kL) 2 ). This factor L 2 is combined 
with the prefactor k^bA 2 in Eq. (|B5j) . From the defini- 
tion of the charge parameter q, we can write the result 
as an overall prefactor qn 2 L 2 . The truncated expansion 
is defined as the expansion up to fourth order in kL of 
the expression where this prefactor is taken out. This 
means that we determine the series expansions of the ex- 
pressions in Eqs. (|D3|) and (|D27|1 . after we divide by a 
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The next domain is Ls ™7 < r < L sin ^ , where the 
expression gets a lot more involved 



K 2 A),o(r;7) = - — 
sin 7 

' [ a(Kr) a (n \ ■ < \r ( rLs ^ 
/ dip Co(«r) - io Kr /Co — — -, — -^y 



kL sin 7 



/■/3(«r) 

/ dipk { K r)l I - 



+ / dtp I C («r) - i (Kr)lC , . 

J/3( K r) \ \^2sm^+^j 

1 cxp(— «r) 



kL sin 7 



4 / 7T 

I 2 arcsm(C) - - 

sin 7 V 2 



sin 7 



(2arcsin(£) - |) 



/■cr 

7r\ sinh(w) 



sinh(Kr) n 9 
H i — - k L 



nr 



7 7T — 7 

In tan — tan 



2 kL sin 2 7 k 3 L 3 sin 4 7 

x I 1 H - 

1 kL 24 2560 



+ ^/l + sin 7(2 - sin 7) — 
+ \/l + sin 7 (16 — 8 sin 7 + 2 sin 2 7 — 3 sin 3 7) 
1 — (7 + 5 cos 2 7) 



^i 3 
3840 



36 



21600 



,3^3 



|arctanh( v ^ j l- + ^ y+ r ;|()() 



2r 



6 1 ^ ' 480 

2 2 



W2f + l)^-(8f« + 4f + 3)^Y 



cxp(-Kr) ^ 



2r 



x ( l + ( 2 ^ + l)— + (8£ 4 + 4^ + 3) 



4 4 \ 
5400 J 



• (E2) 



domain is given by L sin ? < r < L cos ? 



« 2 -4o,o('";7) 



sin 7 



sin 7 



/ dip C (nr) - i (Kr)/C — - 
7o \ \ v 2sm(^+^) 

— / 

f 2 1 , / \^t- / KLsin7 
+ / d(^fco(Kr)X 



2 sin (<p+i) 



4 / . ... 7\ / 1 cxp(— nr) sinh( nr) 

— arcsin(t) — — 

sin 7 V 2/ \nr nr nr 

sinh(«;r) 9t9 
— - k L 

nr 

, / 7\ / 2 «Lsin 2 7 K 3 L 3 sin 4 ' 

- In tan 1 H 

V 4/ V«:i 24 2560 

/I+COS7VKL /I+COS7V. „ . K 3 L 3 

+ (— T^J T + ( 7 - 3cos ^W 

I + COS7 /l + cos7\ 2 _ .k 2 L 2 
' 1 (2 -cos 7) 



1 + cos 7 



36 

« 4 L 4 



2 ; (7-6cos 7 + 2cos 2 7) g4()() 

3„3 



- I arctanh (v/^ 2 ") f - + £ 2 -f + £ 4 * "'' 



nr 



6 



160 



«r ,„ ~ „, « 3 r 3 



- 1 - (2e 2 + 1) 



2 2 

~3lT 



(8^ + 4^ 2 + 3)— 
5400 / 



+ 



exp(— nr) 



k L 



nr 



1 — cos 7 /l — COS7 \ . . /t 

+ I ~ I (2 + cos 7 ) 



2 L 2 



, 1 — cos 7 \ 
+ ( -I (7 + 6cos7 + 2cos 2 7) 



36 

« 4 L 4 
5400 



1 + (2£ 2 + 1) 



K 2 r 2 
~3lT 



(E4) 



We have abbreviated 



Finally, the domain where r > L cos ^ yields a more 
friendly expression 



Lsin7 
2r 



(E3) 



This domain corresponds to the case where the circle of 
radius r intersects the edge of the parallelogram twice 
at each quadrant. The following domain corresponds to 
the case where there is just one intersection per quad- 
rant. Recall that we assume < 7 < ir/2, such that this 



K 2 A,o(n7) 



sin 7 Jo 



d(pko(Kr)Xo 



kL sin 7 
2sin(^+i) 



exp(-Kr) 22 

K L 



nr 



k 2 L 2 / 9 \ k^L^ 

1 + l^ + ( 7 + 5cOS ^21600 



(E5) 



In the case of parallel rods, we can apply the alter- 
native series expansions, or apply the limit 7 — > on 
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the last two expressions above. Both yield the following 
approximations, where for r < L 



/cMo,o(r;7 = 0) 

^ L — r / 1 exp(— Kr) 
r \Kr Kr 

exp(— nr) ( _ 1 , f N 1 cosh(Kr) 



2kL 
2kL 



Kr 
sinh(Kr) 



' / 1 

Kr 



shi(Kr) + 
r(0, Kr) - 



exp(— Kr) 



-T(0,kL) 
exp(— Kr) 



exp(— Kr) Kr / K 2 r 2 
2kL — ^ ^— 1 



2kL 



Kr 
sinh(Kr) 



36 



t<r 

exp(— kL) 
kL 



K 4 r 4 \ 
1800/ 



Kr 



r j kL r 



Kr / Kr 
kL / kL 



36 

k 2 L 2 
36 



K 3 r 3 


K 4 r 4 \ 


240 + 


1800 y 


K 3 L 3 


K 4 i 4 


240 


+ 1800 



(E6) 



For r > L we obtain 



K 2 A,o(r;7 = o) 

expf— Kr) / 1 

= 2kL — — shif/tL) + — 

Kr \ k£ 

expf— Kr) kL / k 2 L 2 



cosh(K-L) \ 
kL / 
k 4 L 4 \ 
1800 J 



(E7) 



Likewise, there are expressions for I = 2,4, These are 
all used together to create an (approximate) expression 
for the pair interaction outside of the hard-core exclusion 
region. We use this pair interaction to numerically calcu- 
late the effective excluded volume. This is accomplished 
by a numerical integration scheme over all different do- 
mains of r, for given rod orientations. Our approach 
is fundamentally different from other theoretical work 
[H, HH , in the sense that we apply the interchange of the 
two positional vectors r and lw — l'ui'. We have to do this 
in order to calculate the full integral over r, in contrast 
to the studies in Refs. [U [HI, where only a description 
is given of the pair interaction for rods at large distances. 
Conversely, if one considers non-spherical charge distri- 
butions on spherical particles, this switch is not needed 
when introducing rotational invariants. 
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